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Abstract 

A 1-D linear gyrokinetic code called AWECS is developed to study the kinetic excitation of 
Alfvenic instabilities in a high-/3 tokamak plasma, with j3 being the ratio of thermal to magnetic 
pressure. It is designed to describe physics associated with a broad range of frequencies and 
wavelengths. For example, AWECS is capable of simulating kinetic ballooning modes, Alfvenic 
ion-temperature-gradient-driven modes, as well as Alfven instabilities due to energetic particles. 
In addition, AWECS may be used to study drift-Alfven instabilities in the low-/? regime. Here, the 
layout of the code and the numerical methods used are described. AWECS is benchmarked against 
other codes and a convergence study is carried out. 

1 Introduction 

In a tokamak, nested closed toroidal magnetic surfaces are used to confine a high-temperature plasma 
consisting mainly of ionized Deuterium. Such magnetized plasmas are known to support various kinds 
of magnetohydrodynamic (MHD) shear Alfven waves (SAW), the properties of which are determined 
by the geometry of the magnetic flux surfaces; in particular, the field line curvature and magnetic shear. 
Resonant and non-resonant interactions between SAWs and plasma particles can lead to excitations 
of SAW instabilities. These instabilities may, in turn, affect particle confinement. In order to optimize 
the tokamak geometry and operating conditions for thermonuclear fusion applications, a thorough 
understanding of SAW physics is crucial. For a review of SAW observations and comparison with 
theory see, e.g., Ref. [1]. 

In order to investigate the linear instability of SAWs, a linear gyrokinetic particle-in-cell (PIC) 
simulation code, called AWECS, is developed. The model equations describe the dynamics local to a 
field-aligned flux-tube using the so-called ballooning formalism. The equations are valid for a broad 
range of frequencies and wavelengths, with focus on temperature- and pressure-gradient-driven insta- 
bilities, while ignoring modes driven by the gradient of the parallel plasma current. For instance, 
AWECS allows to study electrostatic and Alfvenic ion-temperature-gradient modes (ESITG and AITG) 
[2, 3], kinetic ballooning modes (KBM) [4], /3-induced Alfven eigenmodes (BAE) [5], toroidicity- induced 
Alfven eigenmodes (TAE) [6], a-induced toroidal Alfven eigenmodes (aTAE) [7], as well as energetic 
particle modes (EPM) [8] . In addition, AWECS may be used to study drift-Alfven instabilities in the low- 
(3 regime [9]. Here, (3 = 2[XqP/Bq is the ratio of thermal to magnetic pressure, and a = — q 2 Rod/3 /dr is 
the normalized pressure gradient, with q being the safety factor (a measure for the field line pitch), Rq 
the major radius of the torus, r the minor radial coordinate, and Bq the field strength at the magnetic 
axis. 

This paper is organized as follows. In Section 2, we describe the physical model and, in Section 3, 
the numerical methods used to solve the equations. In Sections 4-6, AWECS is benchmarked against 
other codes, followed by a convergence study in Section 7. Concluding remarks and discussions are 
given in Section 8. 
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2 Model 



In this section, we describe the physical model used. After providing an overview of assumptions made 
in the derivation, we describe the equilibrium model, followed by the gyrokinetic equation and the 
electromagnetic field equations. Then the equations are normalized and cast into a form suitable for 
numerical solution as an initial-value problem. For convenience, in the first part of this section, all 
time-dependent variables are Laplace-transformed (d t = d/dt — > —iuj). 

2.1 Assumptions and formal ordering 

We employ the linear gyrokinetic field equations derived by Zonca & Chen [10]. A reduction similar to 
that described in Ref. [10] is applied, except that finite-Larmor-radius (FLR) corrections for thermal 
ions are retained in the present paper. The model is valid for Alfvenic instabilities in a wide range of 
frequencies u> and wave numbers k, provided that 

(i) : w < k\\v te , kj_p ce -C 1, (ii) : lo < w ci , (iii) : w < loa] (1) 

where lj cs = e s B/m s is the cyclotron frequency and p cs = v±/u> cs the Larmor radius for particle 
species s (s — c for electrons, s = i for thermal ions, and s = E for energetic ions). Here, k\\ is a 
short-hand notation for the typical value of (B/B) ■ V = d/dl in regions where the field perturbation 
has a significant amplitude. Correspondingly, k± measures the wave number perpendicular to the 
equilibrium magnetic field B. The thermal velocity Vt s is defined as T s = m s v% s . The restrictions in 
Eq. (1) mean that (i) we consider wave dynamics which are adiabatic with respect to electron dynamics 
and neglect electron-FLR effects (formally: m c — > + ), (ii) the magnetic moment is conserved (no 
cyclotron resonances), and (iii) the interaction with fast magnetosonic waves is negligible. Although, 
assumption (i), uj <C k\\v te , is applicable only to passing electrons, kinetic effects associated with 
magnetically trapped electrons are ignored at this stage. Assumption (ii) implies that the only relevant 
"kinetic effect" is the kinetic compression of ions along the magnetic field, so the particle dynamics 
are described in terms of the parallel velocity v\\ and the position I along a field line. 

In this study, effects associated with the anisotropy of the equilibrium particle distribution /o are 
neglected; thus, we let P = P±_ « Pu . Furthermore, since we are dealing with low-frequency waves, 
uj/k -c c, the displacement current in Ampere's law is neglected (V x B m poj). As is typical for 
tokamaks, the plasma is taken to be sufficiently dense to satisfy the quasi- neutrality condition; i.e., 
k± <C 1/Aoei with Ado being the electron Debye length. 

The equations are linearized by separating the distribution J- s (phase-space density) and the elec- 
tromagnetic fields Etot ari d -Btot m to equilibrium and perturbed components, 

F s = f Qs + Sf s , Btot = V x A tot = B + SB, E tot = 5E = -V8<j)-d t 8A; (2) 

where we have assumed that the equilibrium is static (_E t ot = SE). We use the Coulomb gauge 
V ■ A to t = 0. Since we consider micro-scale instabilities, the length scales of equilibrium and perturbed 
quantities arc disparate in the direction perpendicular to the magnetic field; i.e., fcj^ 3> feo_L- Due to 
the stabilizing influence of magnetic tension, perturbations tend to be aligned with the magnetic field, 
so that kj_ 3> fc|| and the perturbations have a nearly flute-like structure. 

In a tokamak, the toroidal component of the magnetic field is much stronger than the poloidal 
component, so that B f=a Bt- In simple toroidal coordinates (r, i?, Q (minor radius, poloidal/azimuthal 
angle, toroidal angle), Bt has the form 

B T = B /R, with R = R/R = 1 + e cos (3) 

The typical value for the inverse aspect ratio for the flux surface under consideration, e = r/Ro, is 
taken to be of the order e ~ 0.1 • ■ -0.3, depending on the proximity to the magnetic axis. The scale 
length of the equilibrium pressure gradient is taken to be of the order L p /Rq ~ 0.1. Using 5 ~ L p /Rq 
as a small expansion parameter, and noting that ko± ~ L~ , we adopt the following formal ordering 
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[10]:* 

w/u A ~ OiS 1 / 2 ); (4a) 
ko\\pd ~ fcoxPci ~ C(<S), fc||p ci ~0(5), fc_Lp ci ~l; (4b) 

2«5 ~ ©(*»/»), A = J^^ 0((5) , ^ = ^^0(^/2) (4c) 

n oi "^AO m i W AO n Oi 

where «ao = Bq/ y/po m i n oi is the Alfven velocity. Furthermore, using the disparity between electron 
and ion mass, m c jm\ ~ 0(<5 3 ), and assuming T c jT\ ~ 0(1), we have 

«i/Wc ~ J0||i/j0||c ~ ^±Pc ~ C((5 3/2 ) -> .7o 1 1 ~ io||e- (5) 

Thus, we may assume that the equilibrium current is carried primarily by electrons. From the above 
ordcrings (4a)— (4c) it follows that 



uj v A 



LO ci U! A0 V ti 

The density and (3 orderings (4c) imply that 

Ti/T E ~ 0{8), k ±PcE ~ OiS- 1 ' 2 ). 
The wave number ordering (4b) implies that [11] 

ik x ■ SA/6A\\ - 0(S). 

Note that the above constitutes a maximal ordering, chosen to cover a broad range of plasma 
parameters; more typical values for the normalized wavenumbcr used in simulations are, for instance, 
k±p c \ ~ 0.1 and k±p cE ~ 1, which corresponds to TI/Te ~ 0(8 2 ). 

2.2 Equilibrium model 

The separation of temporal and spatial scales described by Eq. (4), in particular, to <C u) C i and fc ~ 
k\\ <C A;j_, allow us to average over the rapid gyromotion and apply an eikonal approximation, which 
simplifies the problem significantly. The eikonal approximation is facilitated by the so-called ballooning 
transform [12]. The essential physical effects arising in toroidal geometry are then captured by the 
so-called s-a model [12]. 

For the equilibrium distribution, the separation of temporal and spatial scales implies that, at 
lowest order, /o s may be taken to be independent of the position along a field line I and the gyrophasc 

£ [ii], 

difvs = drfos = 0. (6) 

Since the field lines cover magnetic surfaces ergodically, Eq. (6) implies that the equilibrium particle 
density no s is a function of the minor radius r only. The plasma particle distribution is assumed to be 
an isotropic (Tj_ = Tj|) equilibrium, and we usually decompose it as /o s = no s (r)Fo s (r,£). Thermal 
ions and electrons are taken to obey a Maxwellian velocity distribution 

f 0s = n 0s (r)F Ms (r,£), F Ms = [2irT s (r)]- 3 ' 2 e- s / T ^; (7) 

where 

£ = (vl+ «jj)/2, T s = k B T s /m s . 

Energetic ions, such as a-particles from nuclear fusion reactions and injected beam ions, are generally 
non-Maxwellian and anisotropic. In the present model, isotropicity is assumed for simplicity, while 
allowing for a non-Maxwellian distribution (e.g., slowing-down distribution). 

Due to axisymmctry, the tokamak equilibrium is effectively a 2-D system, with £ being the ig- 
norable coordinate. Taking advantage of the spatial scale separation, an eikonal approximation would 
allow to decouple the remaining two dimensions to create two simpler 1-D problems; however, its direct 



Ref. [10], the ordering of the perpendicular wave number is fcxPci ~ C(<5). Here, larger values are allowed and 
FLR corrections will be included. 



3 



application is infeasible because the magnetic shear severely distorts the flux tubes, while all physical 
quantities must be 27r-periodic in •& and £. The periodicity constraint is eliminated by the ballooning 
transform, which maps a flux tube of the periodic physical system onto a non-periodic covering space, 
with the new coordinate being the ballooning angle 9 £ (—00,00). The physical solution is recon- 
structed from the solutions on the infinite domain by linear superposition. Formally, the ballooning 
transform in an axisymmetric configuration can be written as [12] 



00 

(r,0) = — e ~ ml? / Me~ ie(nq ~ m) S<i>(r,0) = <^M + 27rp); (8) 

m— — 00 



where m is the poloidal and n the toroidal mode number. The transform is valid for short-wavelength 
modes {nq 3> 1). The dependence on r is then separated as 



6<j)(r,9) = 5§{9)W{r)e lS , 

where S is the cikonal (rapidly varying with r) and W the slowly varying radial envelope. 

In order to proceed further, one has to construct a local model for the MHD equilibrium. A 
popular choice is the s-a model by Connor, Hastie and Taylor (CHT) [12], where we may write 
8(/>(r,6) — 6&(6)e lkrr . In the simple case of a cylindrical tokamak, the radial wave number is k r = 
—ktfs(9 — 9k), where s = (r/q)dq/dr is the global magnetic shear, 9k oc ko± describes the slowly-varying 
radial WKB envelope, and k$ = nq/r is the poloidal wave number. In toroidal geometry, the thermal 
pressure modifies the flux surfaces and, thus, imposes a poloidal modulation on the magnetic shear. 
The CHT s-a equilibrium model retains the lowest-order effect by assuming that the toroidal magnetic 
flux surfaces maintain a circular poloidal cross-section while undergoing an outward shift (Shafranov 
shift). In this case, the gradient and Laplacian operators applied to perturbed fields, and the magnetic 
curvature drift, have the following form: 



fcij ktf 



-9 k )-asm9 = h{9), (9a) 
+ h 2 (9) = f(9); (9b) 



fi K = k x ■ (b x k) = — § -g{9), g = cos9 + h(9)sm9, (9c) 
K{V) 

so that kx = k$\fj and the local magnetic shear is dgh = s — acos(9). Here, b = B/B and k = b - V6 
is the field line curvature vector. In this model, the equilibrium is completely described in terms of 
two parameters: the flux-surface-averaged magnetic shear s and the normalized pressure gradient a. 

Since the CHT s-a model is derived for the low-/?, large-aspect-ratio limit, 2 its application in a 
study of high-/3 instabilities in tokamak plasmas with e > 0.1 is somewhat ambiguous. To justify 
its use, it is argued that it does capture essential effects of toroidal geometry and high /?; namely, 
toroidal curvature and modulation of the magnetic shear. The presence of these features physically 
distinguishes a toroidal plasma from a cylindrical or slab geometry. On the other hand, higher-order 
toroidal effects ignored by the s-a model (such as elliptic and triangular deformation of the flux 
surfaces) merely modify the lowest-order results; for instance, by adding further frequency gaps and 
corresponding Alfven eigenmodes. Thus, the s-a model is thought to be a convenient and powerful 
tool for basic studies of the qualitative features of low- and high-/? tokamak instabilities. 

2.3 Gyrokinetic equation 

The evolution of the particle distribution T s is chosen to be governed by the Vlasov equation; i.e., the 
ensemble-averaged phase-space continuity equation in the absence of collisions. For a high-temperature 
tokamak plasma, where electromagnetic forces dominate, the Vlasov equation reads 



+ v . v + -2- (E + v x B • — 

at m, ov 



Ts = 0. (10) 



2 For a detailed derivation of the CHT s-a model see, e.g., Ref. [30]. 
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Transformation into guiding center coordinates, linearization, and application of the approximations 
and orderings outlined in Section 2.1 yields the collisionless linear gyrokinetic equation (GKE) as 
derived by Antonsen & Lane [11]. Here, we adopt the notation used by Chen & Hasegawa (CH91) 
[13] and neglect anisotropic contributions due to d^fos, where H = v]_/(2B) is the magnetic moment. 
The fast gyromotion and the electron response to the parallel electric field 8E\\ = di(Stp — 6(f)) are 
eliminated through the substitution 



5f„ = — Uudsfo.) W - JlH) ~ (^ s fos) J%H] + JoSK 
m 



-CH91. 



(11) 



where 



£ = v 2 ± /2 + vf { /2, 

didifj = iuSAn (with Coulomb gauge), 
u* s = ^/(fex x b) ■ V g , 

and V g is the Laplacian in guiding center coordinates. Defining 5G S = cuSK^ 91 , Su = u)(8(f> — Sip), 
and applying the gyroaverage, one obtains 



[v\\di - i(u - u) ds )] SG S = %—— (Qfos) {SS 1s - ia5S 2s ) 



with the source terms 



5Su = J 0s 6u + u ds J 0s 5ip + —JisuSBn 



SS 2 



hlKcWi^; 



(12) 



(13) 



where Ji(X s ) is the Bessel function of i-th order, A s = k±v^/uj cs , and a = v\\/\v\\\. The magnetic drift 
frequency is 

LUds = fi re i|/w C s + £Ib^B/lu cs = Q K (vj + fiB)/uj cs + Qpfj,B/uj cs , 

with 

fl K = k± ■ (b x k), n B = fc ± • (b x VlnS) 
tip = -(Mo/B 2 )fc_L ■ (b x VP). 
The phase-space-gradient operator Q = uj&£ + cD* s may be written as 

Q^Q = {u£ s -u® s )/T s , 



where 



Here, 



-u* s L ns d r ln/o s — 



T~2 



T-l _ "OS 

j- ) -^ns i 'is if i 

WcOs-Lns n-Os n 0s/ n 0s 

and the prime denotes a radial derivative. Note that we have used the sign convention £ • b < 0, so 
that 

(fej_ x b) ■ V g « —k®d r . 

Finally, application of the ballooning transform (qR$di — — > dg) and the CHT s-a model yields 

. e s /os 



e. 



Vs = 



-T s d £ Infos ^ 1. 



(14) 



V/T 



—=-dg - i(ui - LO As ) 
qR 



6G» = 



(w„ s - u)Q s ) (5S ls - ia5S 2s )] 



where the coefficients in the magnetic drift frequency uj^s ar e given by 



fi„(0) 



R(6) 



9(0), 



ktfCt 



2q 2 R(d)B 2 {8) 



(15) 



(16) 



For passing electrons, Eq. (15) is dominated by the v^dg term alone, which allows us to set <5G e , P ass = 
0. As mentioned above, in the present work, we ignore trapped-electron effects and set <5G e ,trap = 0. 
Thus, Eq. (15) is solved only for ion species, while electrons are treated as a massless fluid. 
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2.4 Electromagnetic field equations 

The evolution of the electromagnetic fields is governed by Maxwell's equations. We adopt the field 
equations as derived by Zonca & Chen [10], including terms containing deX s and neglecting anisotropic 
contributions due to After application of the ballooning transform, the equations read 

k 2 d ( dSib\ t-^ I e 2 \ 
° = J^Rofde V^f) "^Ey 1 - Jo)Qfo) a "fy - V P [{Sip + 2Sl K )5i> + 

- ^ E (S Wd(1 " J ^ Qfo ) S ^-^Y1 ( e ^ B ( x - J °) Qf°) " 5B \\ 

-f, J2(^dJoSG) s + in J2(e^(d 6 X)J 1 6G'j , (17a) 
0=J2(eJoSG) s + J2(^d£fo) wW-^) + i;(^(l-J?)0/o) <*V>, (17b) 
=w5B|| + Sl p 8i}> + ^ ^e/*B (l - ^J ) Q/o^) S^j + ^Y,\ m ^ B ^Y SG ) > ( 17c ) 



= (^ s - w e s )/r s , (...),= / d 3 v 



where 



Equation (17a) is the so-called vorticity equation and is obtained by substituting the parallel Ampere's 
law into the continuity equation, which, in turn, is obtained by taking the zeroth-order velocity moment 
of the gyrokinctic equation (15). Equation (17b) is the quasi-neutrality condition, which originally read 

s 

and Eq. (17c) is the perpendicular Ampere's law. 



2.5 Definitions and normalization 

The following dimcnsionlcss parameters are used in the present work: 
r L ns d In r d In r 



£ ~ i? ' £ns ~ Ro ~ £ d\nn Qs ' £ps ~ £ 'd\nP s ' 

dlnT s t m sT s Z s n Qa e s 

Vs = T"j , T ss> = 7r~i T ss' = ~ ) Z s = , (18) 

dlnn 0s m s 'T s < Z s >n 0s , e e 

M s = —, k s = ^-^, b s = fk 2 s , \ s = ^f^±. 

Since we are dealing with a quasi- neutral Deuterium plasma (Z\ = 1) containing only a sparse popu- 
lation of energetic particles, we have e n \ — £nc- Note that e„ s = e ps (l + rj s ). 
The equations are normalized in two steps. In the first step, we let 

v „ _ u A _ ""ao^k _ lk 0s g{9) 



<2> = , il KS 

v aqSIp _ o^os / SG S \ _ v' A0 / dG s 



VAO ' WAO ' ^AO^cs T, 1/2 



WAO^cs 2qTs' B 2 {6) I ^Ma J "Os I *Ms 

6* e = ^, 8U C = A^ — , SC.- 6c 



where 



l c m c l c ujAo m c l c WAo^ci 

v AQ = v A /B 1 u AO = VAo/(qRo), B = B/B = R Q /R. 
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Note that 2% = /3{ = 2^P\j and k s = ko s /B. In the second step, we let 

A 2 



^ks I -nvis 



so that the normalized drift frequencies become 



qko s T s 



1/2 



iu ds = n KS (vf + ijlB) + Vl ps ^lB, (21) 



ak 0s T s 



-1/2 

qk 0s f s 1/2 g(0), n ps = -'- 



Furthermore, we have 

A s = \Zb~ s v_L, 9 S = -%ln/ 0s , wf s (F M ) = 2* g [l + 77 s (£ - 3/2)]. 

In the following, for simplicity, we will omit the tildes and hats indicating normalized quantities, except 
for B and k s . 

The second normalization is motivated as follows. First, the values of input parameters determining 
the velocity space domain size to be sampled numerically become independent of the species' temper- 
ature and may, thus, be kept constant once a set of suitable values has been determined. Second, we 
prevent fi re ; and f2 p i from diverging, and F^i from collapsing into a Dirac 5 distribution in the limit 

1 /2 

Note that, in AWECS, % must remain finite; the case % — can only be realized if T { 1 is included 
in the time normalization, which is only done in an electrostatic version of the code (see Section 4). 

A ■ ■ ... —1/2 

The wave number k s is taken to be parametrically independent of T s , which implies that k$ oc T s . 

All moments in Eq. (17) involving the Maxwcllian distribution, (...Fm s ), can be evaluated analyti- 
cally, which gives rise to the following quantities: 

r 0s = e- b 'I {b.), A 0s = l-^, A ls = 1 + - lV (22) 



4 2 \2 Iqs J 2 \Iq s / & \ ios / \ J 0s 

T 0s = A 0s + Vs (1 - 26 s A 0s ) = A 0s + 2(T ls - 1) + 77,, T ls = 1 + r) s b s ( -il - 1 ) , 



T 2ps = ( A is - - ) 



with Ii being the modified Bcssel function i-th order. These quantities are used in the next section. 
2.6 Final form of the equations 

In this section, we write down the normalized equations in a form suitable for numerical solution as 
an initial value problem. In the following, we denote dgf by /' and the Laplace transform in time is 
undone. For energetic ions, the normalized velocity space integrals are given in both the general form 
valid for any equilibrium distribution, Fq s = fos/no s , and the form for a Mawellian (Fq s = Fm s )- 

Marker motion and gyrokinetic equation 

The evolution of the position 9j of a particle labeled j is governed by the equation 

d t ^=T// 2 % .(^). (23) 
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In order to eliminate ui idt on the right-hand side of the GKE, the slow response of the ions is split 
as follows: 3 



SG S = Sg s + Q s F 0s (SSis - iaSS 2s ) ■ 



which yields (for s^c) 



£o8g s = -icodsSgs + iF as [{ujJ s - @ s u<i s )(5Si s - ia5S 2s ) + S s <5Ai s + 6 S 6A 



2s 



Here, Cq — dt + T s 1 ^Vnde is the propagator and the source terms are 

T 2 

6S U = J 0s 5U s + (j ds J 0s 6V s + ^-^Ji s 5C s , 6S 2s = -T s 1 / 2 \v ]l \X' s Ji s S* s . (26) 
The derivatives <5Ai s = iv\\dg5Si s and SA 2s = \v\\ \d$SS 2s are given by 



(24) 



(25) 



SA ls =iT s r /\ I J 5U' + u d J 5*' + u> d J 6* + ^^J X 5C + ^^f^J x 8C 



B' 

6A 2s =T sf iB—\' s J ls 5V s - T s vj I A'Jii*' + A 



T 



A' x 2 



Ji<5* + (A') 2 JoSV 



(27a) 



(27b) 



where SS S = —Z s T^ s SS e (similarly for 6^ s , SU S , and 5C S ), J, = Jj(A), and A s = y/Jk s v_\ 
derivatives of X s and Wds with respect to 9 are 

A' s /A = h(s - acostf)// - B'/(2B), 

K/X s - (A'/A) 2 = [ha sin0 + (2 - f)(h') 2 }/ f + B"/(2B) - (B>/B) 2 /2, 
uj' ds = (Cl KS /g)[hcos6+(ti - 1) sin#](u 2 + fj,B) - (fi res + n ps )nB(B' /B). 



The 



(28) 



Vorticity equation 

The vorticity equation is a second-order differential equation in time t. Collecting all time derivatives 
in an auxiliary field 8E C , we obtain two first-order equations: the continuity equation 4 

d t SE e =k 2 0i + 2hti6& e ] - [u;*; (1 - ToiTn) + H EU ] SU C (29) 

+ w»i [2tJ(1 + 7 h )Q Ki + t£(1 + r) e )tl pi + 2fi Ki r 0i T 2Ki + 2fi pi r oi T 2pi ] <5* c 
+ H E *6* C + w*; [rj(l + Tfe) + r oi Toi] <5C* C + ir £C <5C c 



1/2 



'T- 1 / 2 -n 
2 E r Ei 



+ — (w d Jo<5G)i + -f (w d J £G) E - i-i^- («||A'JitfG) - i E T Ei (vyA'Ji^G) 



and the parallel Ampere's law (which now defines JE'e) 

A u d t 6* e =% [2(n Ki + n pi )(i - roiAxi) - - r 0i ) - - r ffi Tii)] (30) 

+ ii?**(5* c + <5£ c + iAuSUe + i(l - r iAoi)<5Ce + iHy C $C c ; 



3 Before normalization, this transformation reads SG 3 = Sg s — ^f- (dsfos) (55i s — icrSS2a)- 

4 Note that we have combined the drift-kinetic and FLR components of the ballooning terms (oc <5\I/ e , 8C e ) which are 



separate in Eq. (17a). 
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where A u = (1 — Toi) + H u , and the energetic particle terms are 



H E U =^ETEi r iE ((1 _ ^o) w ^o) E ^ETEi r iE w *E (1 - ToeTie) , 

He^St =Z^T^ i T^ 1 ( JgWd^f F ) E 2Tg;tJ»E [^KiToE^ 2/^E + ^piroiT 2p E] , 

,2 Jl J T \ F M „ "T 

w * ^0 ) — *• TEi tJ *EJ- OE J-OE, 



H EC =t^ ( nB 



H^xn =Z E T Ei T? E ((1 — Jo)wd@s-Fo) E — Heu, 



Fm 



2T£(n ri + f2 pl )(l - r 0E A 1E ) - rSiO pi (l - r 0E ) - ff ££ /(F = F M ), 



H*c =T£i 1- 



2 Ji Jo 
A 



T Ei(l ~ ToeAqe), 



^ =^r£rf E ((1 - J 2 )e s F ) E ^ Z E r^r^(l - r 0E ). 
The moments of the transformation 5G — » (5g are: 



( Wd Jo(<5C?-(5.g)) s 



(31) 



/ ITT 



F M 



r 0s [2 + n pi ) a 1s - n pi ] <5c/ e 



(o Ki + n^) As. + 2n Ki {n Ki + n pi ) a 1s - - + -n z Ki 



i 



z, 



T 0s [Sl Ki (l + 2A 0s - 2&iA 0s ) + fi pi (l + A 0s - 26;A 0s )] SC , 



A, 



(32) 



M, 



A' 



A. 



2rf ^ Mr ,Ao.<y*e- 



(33) 



In the low-/? limit we may set 



5Bn = n„ = 0. 



By neglecting the energetic ion terms we then recover the model used by Zhao & Chen, 2002 [25]. In 
this case, the term on the second line of Eq. (32) reduces to 

Q 2 K [A 3 + 2(Ai - 1/2) + 3/2] = 2fi*A a . 

Note that our calculation yields a different A 2 than that given in Ref. [25]: A| hao = 5/2 + whereas 
here A 2 = 7/4 + ... [cf. Eq. (22)]. Our A 2 gives slightly larger growth rates, as will be shown in the 
benchmark in Fig. 6 below. 

Quasi-neutrality condition 

After eliminating dtS^ c with the help of Eq. (30), the quasi-neutrality condition becomes an algebraic 
equation, which may be written as 



1 



SU C = - iSE e + [2(fi Ki + f2 pi )(l - ToiAn) - n pi (i - T 0i )] 6V C 
+ H U9 S^f e + (1 - T 0i Aoi)SC c + H UC SC C 



(JoSG), - -f (JqSG), 



(34) 



The energetic particle terms are 



H V v =Z E T£ i T? E ((1 - J 2 )^ d eF ) E ^ rg; [2(fi rei + a pi )(l - r 0E A 1E ) - M 1 - Toe)} 
2 Ji Jo 



Hue =r£ 1 



A 



LiBOFo) ^ 7^(1 -IWW 

/ E 



(35) 
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The moments of the transformation SG — > Sg are: 

(J (6G-6g)) s Z s tI { , 2 ^> :tt , . „ . , 1 / 2J lt / ( 



[< J 2 9F ) s fl7 e + (w d J 2 6F > a 5V e ] + ( M B^-^9F ) SC, 



^ Z s r^r 0s( 5(7 c + T 0s [2tt Ki A ls + (2A ls - 1)] <5* c + r 0s A Qs SC c . (36) 

Given the orderings described in Section 2.1; in particular, ~ ~ 0(<5 2 ...<5), we expect the 
contribution of energetic ions to the quasi-neutrality condition to be small. 



Magnetic compression 

After eliminating dt5i&c with the help of Eq. (30), the perpendicular Ampere's law becomes an algebraic 
equation, which may be written as 

(b 2 + %AlAS) SC c = - i%A^5E c + TiAxAuSUe + 7IA E iW#e 

+ %a^ [2(n Ki + n pi ){i - roiAii) - n pi (i - r oi ) - ^(1 - r iT u )] 

+ Tiu^(l + 7] C )S^ C + TiUj^Toi [A oi + rn (1 - 26iA oi )] <5* c + i?c*<5*c 
+ — ^ (v^SG), + _^ T \ r <ttx-W E ; (37) 

where 

j4 s = -[(i-r Qi Aoi) + J?£tt 7 ]/A,. 

The energetic particle term i?c* is 



Hc9=-%Tg l ^—^ l iBu,iFoJ ^ 7Ir^ E r 0E [A 0E + ?7e (1 - 2fe E A 0E )] ■ (38) 
The moments of the transformation 5G — > (5g are: 



rr T I oj2 



v5J Z s \ A 
^ - ^r 0s [fi rei (l + 2A 0s - 26 s A 0s ) + O pi (l + A 0s - 26 s A 0s )] <5*c 
27V T 

- -^r 0s A 0s 5C o - 7ir 0s A 0s( 5C/ c . (39) 
The final form used in AWECS is obtained by substituting Eq. (37) into Eq. (34) to eliminate SC e . 



3 Numerical methods 

Equations (25), (29), (30), (34) and (37) are solved as an initial value problem with the particle-in-cell 
(PIC) code AWECS using a Runge-Kutta scheme. A finite number of markers is employed to sample 
the phase space. A modified Sf method appropriate for particle-conserving compressible dynamics is 
adopted, with marker weights chosen such as to allow a uniform distribution in energy. In this section, 
these methods are described in detail and an outline of the computational cycle is given. In the 
following, grid points and markers arc labeled by the indices i € [1, N g ] and j g [l,JV m ], respectively. 
A "cell" is the space between two grid points and its size is A9 = — 0i. 



3.1 PIC method 

While marker positions 0j vary continuously, fields are sampled at discrete grid points 9i. In order to 
solve the field equation, the contribution of each marker to the particle density at each grid point must 
be determined. Conversely, the marker motion is subject to electromagnetic forces known only on the 
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discrete grid. The PIC method employed here is a lst-order scheme that provides smooth mapping 
between markers and the grid. In this method, each marker j is replaced by a top-hat function II of 
width AO centered at Oj ,where the definition of the top-hat function is H(x) = 1 for \x\ < 1, and zero 
elsewhere. The sum of these finite-sized markers integrated over the interval £ [0» — A9/2, 9i + Ad/2] 
and divided by AO yields the number of markers Ni contributing to grid point 9%. 



A6 



N 6,+A6/2 



J ~ ± di-Ae/2 

This local integration gives rise to the triangular shape function 5 

S (9-0.) = ! l-W-W/M ■ \0-9 j \<M/2, 
' 3 ' I : elsewhere. 

Note that J d9 S = A9. At the boundary points of the domain in which markers are loaded, the 
marker weights are doubled. This effectively simulates the effect of a plasma beyond these points, 
which is a mirror image of the plasma inside. 

3.2 Modified Sf method for compressible dynamics 
Description of the method 

Equation (15) may be compactly written as 

CoSGs + iwdsSGs = SCfos, 

where Cq = d t + V\\di is the inverse propagator along a field line. The lj^s term can be eliminated by 
a transformation from the guiding centers to magnetic drift centers [4] 



SG S = 5G ds e- i5d ° with <5 ds = / dr w ds (r) 



C 5G ds = -6Cf 0s e iS *°. (42) 

The integrating factor exp(— iSds) was expected to help avoid numerical problems associated with the 
secular term in uj^s {u& s 9 for \9\ 3> 1). To date, however, no significant difference was found 
between AWECS runs solving Eq. (15) and those solving Eq. (42). 

For an isotropic Maxwcllian equilibrium distribution, /o = no{r)F^{r, £), the particle density at 
high energies £ is low, so a corresponding marker distribution would introduce a large amount of 
discretization noise through the highest-order velocity moment, which in our case is (u>d.Jo5G). One 
way to avoid this problem is to load the markers with a probability distribution function (PDF) 
defined by PaWo = /o, and require that dsPa = 0. For the perturbed particle distribution this means 
that PaSW = (5/, or, equivalently, 5W/Wo = Sf/fo, where W and SW are equilibrium and perturbed 
weight functions. 

In the conventional Sf method, the PDF P& for the markers is chosen to be such that L^P„ = 0, 
and one solves the equation for the weight function SW instead of Sf. In the present case, that equation 
would read 

[Co - {div\\)]5W = -W SC(\nf 0s ). 

In AWECS, we adopt the alternative scheme utilized previously in Ref. [14], where P„ is required to 
satisfy the continuity equation for compressible dynamics, 

d t P a + dg^Pa) = 0. (43) 



5 It may be more intuitive to refer to II as the "shape" of a particle. The function S describes the mapping between 
marker position and the finite-difference grid and is sometimes called "assignment function." In practice, we only deal 
with S and it has become customary to simply call it "shape function." 
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Defined this way, the spatial marker distribution represents closely the physical particle distribution 
along a field line. Therefore, we can solve the GKE for the physical particle distribution function, 

£ 5G d s = -SCf 0s exp(i6 ds ), 

along unperturbed marker orbits. The difference between markers and physical particles manifests 
itself only in the discretized velocity space integrals, where an additional weight factor appears (see 
Section 3.2). 



Marker loading 

It is convenient to introduce the pitch angle variable A defined as 

A = = V ' 2± - — = sin2 ^) - (44) 
£ v 2 B B ' 

where tp = (f(B) is the pitch angle; i.e., the angle between the local magnetic field B and the velocity 
vector v. Note that for particles which are trapped in a magnetic mirror, the pitch angle coordinate 
A is related to the turning points ±6^ ("bounce angles") through A = 1/B(9^) which follows from 
u l($b) = 2£. For the purpose of marker loading, we take A to be an independent velocity space 
coordinate instead of the magnetic moment /z. Recall that 

£^v 2 /2 = (vl+vj)/2, (x = v 2 ± /(2B), 

so the velocity space coordinates do not contain the particle mass. For the parallel velocity f|| of 
a marker with (A,£) at a given location 6j, parametrized by B = B(6j), we obtain the following 
expression: 



— v\ = \p2£\j 1 - AB = \/2£u, with u = \J 1 - AB. (45) 

Note that the normalized parallel velocity u satisfies dgu = 0. 

In order for the plasma conditions to remain constant in time, the markers must be loaded in 
equilibrium: dt Pa = 0. Thus, Eq. (43) implies di(v\\P„) = 0; that is, 

v {l P a =C(A,£). 

Since dgA = 0, a simple PDF which satisfies the condition vnP& = C(A,£) and allows to initialize 
uniformly in energy {psPa = 0) is obtained with the choice C(A,£ ) = CqV2£: 



P,(A,B) = ^ = ^£^ = ^. (46) 
hi I Vl-AB u 

The size of the computational domain, 9 G [— 6* max , max ] must be sufficiently large to avoid unphys- 
ical reflections at the boundaries. However, markers are only required in the region where the unstable 
modes have a significant amplitude. In AWECS, the parameter N p determines the number of loading 
periods, and the size of this domain, 2irN p , can usually be chosen smaller than the computational 
domain, 2# max (see the convergence study in Section 7). In the region without markers, we set SG = 0. 

The markers are loaded according to the following procedure: 

1. Spatial loading. Integrating the marker PDF Pa- given by Eq. (46), and using J2a ^ = 2-P<r, we 
obtain 



E 



A max e 2 2 

2 



J d£ J &A j d9Pa = 2C Q (£ max £min ) / d0 



1 - AB 

B 



(47) 



which is valid for any interval [9±, O?]. Hence, the marker density n(9) is given by 



n(9) oc w{9) = I (J 1 - A~B~ - ^max{l - A max B, 0}j = Mmax ^" min . (48) 

The weight function w(9) is then used to map a uniform distribution of random numbers Rj € 
[0, 1] to a nonuniform distribution in 9: 
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(i) Numerical integration of w gives the cumulative distribution 

W{9) = / d8'w{6'). 

Je-L 

(ii) From this, we sample uniformly distributed random values, 

R j W(6 2 )=W j e [o,W(6 2 )}. 

In ballooning space, the integral limits are 8\ = and 82 = max{#b} — Enum for trapped 
particles, while for passing particles 82 = 7r — £ num . Here, e num is the smallest number that 
can be represented numerically in double precision. 

(iii) The map W {8) : Wj -> Oj then yields non-uniformly distributed marker positions 9j. 

Since W is known only at discrete grid points, the inversion W is done using linear 
interpolation. 

(iv) Offsets ripj-vr, with n P j = 1, (N p — 1), are added to spread markers over all periods. 

2. Pitch angle distribution: The largest allowed value for A depends on the location 8 (parametrized 
by B). The constraint to be obeyed is most easily written for the parallel velocity, which must 
satisfy 

v 2 - v]_ = v\ > 0. 

Thus, we first determine the limits of u = \vu \/\f2£ for a chosen interval L4 m i n , A max ], using the 
relation 

u = \Jl- AB. 

These are 

Umax = \A - ^min-B, u min = sj max{l - A max B, 0}. 

Next, we sample random values Uj, uniformly distributed over the interval [u m in, it max ], and 
calculate the pitch angle variable Aj = (1—Uj)/B. Note that passing particles satisfy A max = A n , 
where A„ = 1/B(n) = 1 — e. 

3. Energy distribution: We distribute marker energies uniformly in a chosen interval [£ m in, £max]- 
The appropriate limits depend on the problem at hand; in particular, the order of the highest 
energy moment. In AWECS, the energy coordinate £j is also used to store the sign of the parallel 
velocity, o = sign(uii). 

After loading N ms /4 markers onto the positive 8 and £ axes as described above, the distribution is 
copied to the respective negative axis. The number of markers used for species s is thus given by 
JV ms = Nf s x 4 x 7Vp, with input parameters N{ a and N p . The remaining input parameters specifying 
the marker distribution are 

^mim ^raaxi <C ££71-111111 ^ ^7rmax ^ lj <C $b,min ^ ^b,max ^ ^ : (^^) 

" v ' y ' 

passing trapped 

such that 

G [1/B(8b, max ), l/i?(#b, min)]- 

The fraction of trapped particles can only be manipulated through the inverse aspect ratio e. However, 
the number of markers used to sample the phase spaces of trapped and passing particles can be varied 
independently. 

The marker loading is completed with the calculation of the normalization constant Co- If we 
equate the integral of Pg, in Eq. (47) with that of the Klimontovich distribution, 

iV m 

P s =J2S(£- £j)S(A - Aj)6(8 - 8^)5^ - Zj), (50) 
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we obtain the following equation for the constant Cq\ 



Co = I(f — ~f — \h Wlth bl= / de 5 ; (5l) 

^\^max ^min^l J D 

6»mi„ 

where N Ta (9 m i n ,9 Tnax ) is the number of markers in the interval 9 G [# m in, tax]- In the PIC method, 
the spatial Dirac deltas 8(9 — 9j) in the Klimontovitch distribution (50) are replaced by finite-sized 
markers S(9 — 9j)/A9 [cf. Eq. (40)]. The local normalization constant Coi at a grid point i is obtained 
by carrying out the integral in Eq. (47) over the spatial interval 9 £ [0, — A6/2,6i + AO/ 2] and 
substituting Eq. (40) for the left-hand side: 



^max ^min 



S(9, -9j)=Ni!» 4C 0l (£ ma x - Sn 

3=1 



D 



A9- (52) 



where B was taken to be constant within the small integration interval. As in Eq. (40), Ni denotes 
the effective number of markers at the grid point i and we write 

Cq . = ^lEl . (53) 

^A9(£ max ^min) [^max ^minj.j 

The normalization constant is then obtained by averaging over all populated grid points: 

°o = ^2°oi I ^2^n,>o; where 6 Ni>0 = I Q \ * lge ' (54) 
t=l / i=l L ' ' 

Alternatively, we could define 



C A9 = l j^fmax - fmin) 



^max ^min)/-^]'i J i 

which gives a very similar value. However, the method of calculating local values Coi allows to check 
whether Coi is indeed approximately the same everywhere. 

Velocity space average 

The velocity space average of a quantity X = X(£, (jl, 9) is given by 

oo E / B 

(X)(9) = 2TrJ2j d£ J dfi-^-X^^e). (55) 
a o o " 

Substituting the pitch angle variable A = [iB^/E for /i, and the normalized velocity u = \vu\/"\/2£ for 
and sampling the phase space by including the factor Pg/P^ into the integrand yields 

oo l/B 

(X) (9) «2tt£ Jd£ J dA ^~-X(E 1 A,9) 

a o o " 

oo l/B 

tJ2 f d£ I dABV£PsX(£,A,9); (56) 



2- 



V2C a 

" 

where P& = CoV2£ /\v\\ | was used. Replacing the Dirac distribution 8(9 — 0j) in Eq. (50) by the shape 
function S(9 — 9j)/A9, we obtain 

2 N m 

(X)i « ^ CoAe E B(9 3 )^£~S(9t - Bj)X(£ h A h fl 3 -). (57) 
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3.3 Boundary conditions 



The size of the simulation domain, 2# max , depends largely on the magnetic shear s. The larger s the 
faster any wave with finite radial extent is damped with increasing \9\. However, some eigenmodes, 
such as oTAEs, couple to the Alfven continuum, which consists of waves singular in r. Since r and 
9 arc related though a Fourier transform, continuum waves in 0-space take the form of undamped 
outgoing harmonic waves. In order to prevent unphysical reflections at the boundaries of the finite 
computational domain, the outgoing-wave boundary condition is facilitated by a "boundary filter" of 
the form 



which provides for artificial damping. Typically, 8^ = 0.8 x # max is used. 

3.4 Filtering 

We use the fftw package (version 2.1.5) for optional low-pass filtering in fc-spacc. If enabled, the 
default cut-off is one quarter of the Nyquist wavenumbcr. In principle, this setting may reduce the 
accumulation of aliasing errors, but it cannot eliminate them. In linear gyrokinetic simulations, the 
amount of aliasing is mainly determined by the number of phase space markers and the shape function. 
To date, no significant effect of filtering on AWECS simulation results has been observed, except for 
smoother mode structures. 

3.5 Algorithm 

The time integration can be carried out using either a 2nd- or a 4th-order Runge-Kutta scheme [15]. 
The results presented in following sections were obtained with the 4th-order scheme. The computa- 
tional cycle may be outlined as follows: 

1. Solve the GKE CoSg = —6Cfa s exp(i5d s ) along unperturbed marker orbits. 

2. Push markers along unperturbed particle orbits: 8j(t) = J dt' v\\j(Qj(t')). 

3. Calculate velocity space moments of the marker distribution. 

4. Evolve the electromagnetic fields SE C and S^ e , given the moments of Sg. 

5. Solve the algebraic field equations for 8U e and <5C e , given SE e , S^ e and the moments of Sg. 

AWECS automatically parallelizes on clusters using Message Passing Interface (MPI) . A two-dimensional 
Cartesian processor grid is established in order to allow parallelization over markers and cases. Thus, 
AWECS is capable of parallelizing case scans and producing organized output without the use of separate 
script files. 

4 Benchmark 1: Electrostatic instabilities 

4.1 Preliminaries: ESITG equations 

In order to test the correct implementation of the solver algorithms for field and particle dynamics 
as well as the PIC method, a simple model for linear electrostatic ion-tempcraturc-gradient-driven 
(ESITG) modes is implemented in AWECS-esitg. After renormalizing the time as t = tui*i, with 

1/9 

= gfeoi^T /e n , and writing 




|0| < 0bf, 

|0| > e h{ -, 



(58) 



8<j) 



T?8U e /uj = ei6(f>/(miTi), Sh^Sgi/uj, 
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we obtain the electrostatic limit = SC C = £l p = 0) by letting % — ► 0. Energetic particles are 

excluded and the resulting model equations are 



54> 



l + rT(l-r oi ) 



(JoiShi) , 



dt qk oi 
ddhj 
~dT ' 



itJdiShi + iFi 



Mi 



(cD^i -a}di)^oi<5^ + »^-W||9e(J 
qkoi 



(59) 
(60) 
(61) 



where 



= l + ?7i (£-3/2) , d) di = £„ 5 (vi/2 + , e„ =£ pi (l+77i). 



Note that, in the electrostatic limit, <5/i; = Sgi/u in Eq. (61) is related to the perturbed distribution 
6fi through the relation 

Sfi = ~F m (l - J odH + JoiSh, 

± i 

so it has a clear physical meaning: JoiShi is Sfi minus the density perturbation caused by the ion 
polarization drift. 



4.2 Zero inverse aspect ratio 

In order to test the correct implementation of passing particle dynamics, parameter scans with respect 
to rj[, ktfp s = (r^Y^koi, e n and are carried out for cases with e = 0, studied previously by Dong, 
Horton & Kim, 1992 (DHK92) [16]. DHK92 solve an integro-differential formulation of the model as 
an eigenvalue problem. The resulting growth rates 7 and frequencies u> r are shown in Fig. 1. Zhao, 
2001 (ZhaoOl) [17] repeated DHK92's calculations with a code that solves Eqs. (59)-(61) as an initial 
value problem using a predictor-corrector scheme. Both calculations gave very similar results, so we 
have plotted them as a single curve for each case in Fig. l(a)-(c). A noticeable difference can be seen 
only in the scan shown in Fig. 1(d). The parameters used are 

• Physical parameters: e = a = 0. All other parameters are shown in Fig. 1. 



• Numerical parameters: N m = 2048 x 4 x 3, 6* max = 20, _/V g 



following sections, we use u m i n = 0.01, w max = 5.0, a 
unless otherwise specified. 



256, At = 0.2. In this and the 
in = 0.0001, (Wa* - 0.9999 [cf. Eq. (49)], 



Note that DHK92 defines the thermal velocity as %° ns = V^li so we have scaled their input pa- 
rameters and results by a factor l/v2 where appropriate (fc^pP° ng = s/2k$p s ). Furthermore, we have 
defined the diamagnetic frequency such that w*i > 0. 

Figure 1 shows that the results by DHK92 and ZhaoOl (squares) are reproduced accurately by 
AWECS-ESITG (circles). Runs with 2nd- and 4th-order Runge-Kutta schemes gave identical results. 



4.3 Cyclone base case 

In order to test the correct implementation of trapped particle dynamics, we compare AWECS-esitg 
results for the Cyclone base case with those produced by two other codes: an electrostatic version of 
the gyrokinetic fully- implicit initial- value code GS2 [18], and the global gyrokinetic toroidal particle 
code Gt3d [19]. In the reference cases used, both GS2 and Gt3d were run in the massless-electron 
limit. Both GS2 and AWECS employ the s-a model and the main difference between the codes lies in 
the numerical scheme; thus, the results should be quantitatively comparable. When comparing with 
the global code Gt3d, we only look for qualitative similarity between the results. The parameters are: 

• Physical parameters: The Cyclone base case parameters are e = 0.18, r]i — 3.114, e n — 0.45, 
s = 0.776. ..0.796, q = 1.4, = 1.0 [20]. The actual input parameters used in GS2 and Gt3d 
differ slightly and we use their settings for our calculations. The parameters to be varied are fcoi 
and rji. 
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(b) tv=2.5, £ =0.2, t:' =1.0, s=0.75, q=1.5 
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Figure 1: Benchmark of AWECS results for ESITG modes against Refs. [16, 17] (inverse aspect ratio 
e = 0). Except for the scan (d), the results by Dong et ai, 1992 (DHK92) [16] and Zhao, 2001 
(ZhaoOl) [17] are very similar and are plotted as a single curve (squares). 
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Figure 2: Benchmark of AWECS results for ESITG modes against Cyclone base case calculations 
performed with GS2 [18] and Gt3d [19] (inverse aspect ratio e = 0.18). 



• Numerical parameters: N m = 2048 x 4 x 3, # max = 20, N g = 256, At = 0.2. The bounce angle 
range is 0b,min = 0.1° and b .max = 179.9°. 

Figure 2(a) shows ESITG growth rates 7 and frequencies aj r in dependence of fco;. The results show 
good agreement between GS2 and AWECS. Figure 2(b) shows results of an 77; scan. Both growth rates 
and frequencies agree with GS2 results. The comparison with Gt3d shows surprisingly good quantita- 
tive agreement in the growth rates 7. The systematic discrepancy of about 20% in the frequencies w T 
may be due to the fact that AWECS is a local code and uses the CHT s-a equilibrium, whereas Gt3d 
is a global code. 



5 Benchmark 2: Shear- Alfven instabilities 

5.1 Preliminaries: MHD SAW equation and terminology 

Let us first consider shear Alfven waves with w ~ <^ao, ignore the temperature gradient (r/ s = 0), 
and treat the plasma in the ideal-MHD limit (6U C = 0). In the cold-ion limit (w r ^> {k\\v\\,LO*i,io<n}, 
^oe = and fcoi <C 1), the contributions from kinetic terms and SC e can be neglected and the field 
equations reduce to the low-/? ideal-MHD SAW equation 

dg(fd g S^ c ) + (1 + 2e cos 6) fuj 2 5V c + agS^ c = 0. (62) 

FLR effects may be taken into account by retaining the correction to the inertia term, which yields 
an MHD SAW equation applicable to a warm plasma [21]. In this case, uj 2 in Eq. (62) is replaced 
by (jj(lo — The substitutions S^ s = y/f5^> c and fl = ui — w* P i/2 turn the SAW equation into 

Schrodinger-like form, 

d 2 g 6V s + (1 + 2 £ cos0)(f! 2 - n2)£* B - VWtftt. - 0; (63) 

where 

VbaU = (s -acos0) 2 /f 2 -acos6/f 
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Figure 3: (a): s-a diagram showing (SI) is the first stable domain, (S2) the second stable domain, and 
(MHD-BM) the high-n MHD ballooning unstable domain for 9^ = 0. (b)-(g): Ballooning potential 
Vbaii(0|s, a ) anc ^ m °de structures of instabilities used for benchmarking. See also Table 1. 



a « «cr!t 


<r 1 


(i) < < ii 


■> (ii 

a <, "crit 


ESITG 


AITG 


KBM 


aTAE 



Table 1: Instabilities and MHD eigenmodes observed in the benchmark simulations. ESITG: elec- 
trostatic ion-temperature-gradient driven mode; AITG: Alfvenic ITG mode; KBM: kinetic ballooning 
mode; aTAE: a-induced toroidal Alfven eigenmode. See Fig. 3(b)-(g) for corresponding mode struc- 
tures. 



is the ballooning potential and fi* = w * P i/4 is a constant offset. Equation (63) describes the propaga- 
tion of periodically driven waves (ecosO factor) in the potential Vbaii- Using an MHD shooting code, 
we solve Eq. (63) to obtain MHD SAW frequencies in the high-/? regime, which will be useful to verify 
results of gyrokinetic simulations. 

A derivation similar to that giving Eq. (62), but for marginally stable modes (7 = 0) with frequen- 
cies ui r — u>*i (again letting rfr = 0), yields the Connor-Hastie- Taylor (CHT) MHD ballooning equation 
[12] 

d 2 e 5^ s - y ba ii<5* s = 0. (64) 

Equation (64) determines the stability boundaries of ideal MHD ballooning modes as shown in the s-a 
diagram in Fig 3(a). Ballooning modes (BM) are short-wavelength (high-rt) pressure-gradient-driven 
instabilities similar to the Rayleigh- Taylor interchange instability. They are localized in regions where 
the field line curvature k is unfavorable; i.e., where n has the same sign as the pressure gradient. The 
s-a diagram in Fig 3(a) shows the stability boundaries a^ it (s) and a^(s), which divide the plane 
into the first stable domain (SI), the second stable domain (S2), and the MHD ballooning unstable 
domain (MHD-BM). The diagram shows the case 6^ = 0; the boundaries are modified for finite Ok [22], 
which is not considered here. The Mercier criterion, which determines the minimum values of s and a 
for interchange instability to occur, is not contained in the CHT s-a model, so that the (MHD-BM) 
domain in Fig 3(a) extends to s = a = 0. 

When kinetic particle compression is taken into account, both stable domains (SI) and (S2) become 
populated with temperature-gradient- and pressure-gradient-driven kinetic instabilities. AWECS has 
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shooting code 


ATAE 


AWECS 


"crit • 


0.608 


0.61 


0.61 


(II). 

a crif 


2.605 


2.61 


2.62 



Table 2: Benchmark of MHD ballooning stability boundaries. 



been developed to study these modes. In the following, we will consider several previously studied 
cases to verify that AWECS reproduces earlier results correctly. The instabilities considered are listed 
in Table 1 and examples for the respective mode structures are shown in Fig. 3(b)-(g). 

The instability shown in Fig. 3(g) was discovered by Hirose et al. [23]. The mode structure peaks 
in regions of unfavorable curvature outside the central potential well, so this kind of instability was 
called "higher-order ballooning mode." Based on the fact that the frequency u r approaches a nonzero 
value as k 0i — >■ (cf. Fig. 15 in Ref. [24]), we conclude that these modes are closely related to discrete 
MHD Alfven eigenmodes known to exist in the second MHD ballooning stable domain [7]. These are 
called aTAE, so this name is used in the present work. The names "ballooning mode" and "AITG 
mode" are reserved for instabilities with, in the incompressible thermal ion limit, w r (A;oi — > 0) — > 
residing inside and outside the (MHD-BM) domain, respectively. 

In this section, we benchmark AWECS for shear-Alfven instabilities in the absence of energetic 
particles. Benchmarks including energetic particles are presented Section 6. 



5.2 Ballooning stability boundaries 
Ideal-MHD limit 

In Fig. 4 we compare results results obtained with a shooting code and two initial value codes, AWECS 
(4th-ordcr Rungc-Kutta) and ATAE (2nd-order leap frog). Both ATAE and the shooting code solve 
Eq. (62). The size of the computational domain in the initial value codes is max w 6tt...17tt, whereas 
the shooting code is run with # max up to 2567T. Thus, the latter is expected to give the most accurate 
results for a£) t or a^j] , since the mode structure becomes very broad in 9 when a approaches marginal 
stability. The parameters used in AWECS are: q = 1.5, s = 1.0, e n = 0.2, k oi = 0.001, = 1, r] s = 0. 
The effects of 6B\\ and f2 p are negligible. Kinetic terms are turned off by setting <5G; = 0. 

Results of scans of the parameter a near the first and second ballooning stability boundary are 
shown in Fig. 4, and the values obtained for a^) lt and a^ t are summarized in Table 2. From the 
good agreement between the three codes it can be concluded that AWECS reproduces the ideal-MHD 
stability boundaries correctly. 
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Figure 5: Qualitative behavior of ballooning stability boundaries for finite fcoi with kinetic terms turned 
on (GK) and off (MHD). The growth rate 7 (a) and frequency cj r (b) are plotted as a functions of a 

near a;5t f° r s = 1-0- 



Finite ion Larmor radius and diamagnetic drift frequency 

It can be shown that for 77; = 0, all FLR terms cancel for modes satisfying w m w*i = w» p i, and the 
CHT ballooning equation (64) is valid for any value of fcoi [25]. On the other hand, for the general case 
lo t 7^ w* p i, FLR effects were shown to have a stabilizing effect on ballooning modes, yielding somewhat 
higher values for a£.j t [21]. 

To verify that AWECS reproduces this qualitative behavior, we repeat the calculation of Section 5.2 
with fcoi = 0.001 and fcoi = 0.212, while still setting ij s = 0. Each of these two cases is calculated once 
with the kinetic terms turned off ("MHD") and once turned on ("GK"). The results of an a-scan near 
the first stability boundary o£ri t are shown in Fig. 5. 

In the case (MHD, fcoi = 0.001), the ideal MHD stability boundary is reproduced with good 
accuracy: a£.- t « 0.615. In the case (MHD, fc i = 0.212), a somewhat larger value a£| t « 0.64 is 
obtained, and the mode has a frequency o; r ~ w*i/3 near the stability boundary. 

Let us now consider the two cases where the kinetic compression terms arc retained ("GK"). 
Extrapolation of the averaged growth rates suggests that the modes become stable near the ideal 
MHD value a£. it ~ 0.615. However, the signals are oscillating which gives rise to the large error bars 
shown in Fig. 5(a) and (b). The oscillation amplitude is larger for fcoi = 0.001 than for fcoi = 0.212. 6 
Note that, in contrast to the case (MHD, fcoi = 0.212), the frequency in case (GK, fcoi = 0.212) 
approaches Lu* p i near the stability boundary. However, this should be considered coincidental; the two 
frequencies happen to be the same for this particular value of the wavenumber, fcoi = 0.212, while 
they differ for different values of k 0l . Yet, for this particular case, the extrapolated a£] t coincides with 
the MHD value, so the result is consistent with the requirement that a£j t be independent of fcoi if 
uj — w*i — uJ* P i{i]i = 0). 

It can be concluded that AWECS correctly reproduces the qualitative behavior of the first ballooning 
stability boundary for finite values of the ion Larmor radius. 



6 In Section 5.3, we will show that the growth rates in the latter case, (GK, fco; = 0.212), compare well with the results 
obtained by Zhao & Chen [25]. 
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first stability boundary, k Q .=0.212 




Figure 6: Benchmark in the first ballooning stable domain and near a£ rit . Growth rates (a) and 
frequencies (b) obtained using AWECS are compared with results obtained by Dong et al., 1999 (DCZ99) 
[26] and Zhao & Chen, 2002 (ZC02) [25]. Two cases arc shown: r)i = and rfi = 2.5. 

5.3 Low-/3 instabilities 

The kinetic excitation of electrostatic and Alfvenic instabilities in the first ballooning stable domain 
and near ct£' it has been studied by Dong, Chen & Zonca, 1999 (DCZ99) using an eigenvalue solver 
[26], and by Zhao & Chen, 2002 (ZC02) using an initial value code [25]. In this section, we compare 
AWECS results with data obtained in these two earlier studies. The parameters used are 

• Physical parameters: q = 1.5, s = 1.0, e„ = 0.2, A; oi = 0.212 (fcdp? ong = s/2k oi = 0.3), r£ = 1.0, 
771 = ?7 C , e = 0. 

• Numerical parameters: N m = 512 x 4 x 5, max = 40, iV g = 512, At = 0.02. 

AWECS was run with and without magnetic compression and identical results were obtained, which 
justifies the approximation 8B\\ = Sl p = used in Refs. [26, 25]. 

The results are shown in Fig. 6. The growth rates in both cases, 77; = and 771 = 2.5, agree 
well with ZC02's results 7 and give the same critical a values: a£.j t (?7i = 0) w 0.61 (KBM) and 
a cr!t( 7 ?i = 2-5) ~ 0.5 (AITG). The small discrepancy in the growth rates of the AITG branch is due 
to the fact that our calculation yields a different coefficient A2 than that used by ZC02 [see Footnote 
before Eq. (32)]. If Af C02 is adopted, AWECS reproduces their result exactly as can be seen in Fig. 6. 
The AITG frequencies obtained by ZC02 agree with AWECS results, except for one point: a = 0.75. 
The reason for this deviation is not known. As mentioned in the previous section, the fact that the 
frequency w r approaches near the stability boundary is coincidental and depends on the value of 

fc i- 

The AITG growth rates by DCZ99 differ from AWECS results and those obtained by ZC02. On 
the other hand, the frequencies agree well with those obtained using AWECS. The reason for the 
discrepancy in the growth rates is not known. ZC02 speculated that the problem may be due to the 
stronger sensitivity of the eigenvalue code to boundary effects. However, the mode structure is broad 
only near the stability boundary, while the discrepancy persists for a > a£) t (see also Section 5.4). 
Moreover, Dong et al, 2004 (DCZJ04) [27] revisited this calculation with (9 max w 23 and reproduced 

7 The growth rates in Ref. [25] seem to be incorrect; a factor of V2 had to be applied in order to obtain agreement 
with AWECS results. The authors of Ref. [25] have scaled the results by Dong et al. [26] by \/2 for comparability and may 
have accidentally scaled their own growth rate data as well. It is not clear whether the frequencies are also affected; in 
Fig. 6, the frequency values from Ref. [25] are plotted without additional scaling, because the frequencies of the ESITG 
branch coincide as they are. 
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Figure 7: Benchmark for s = 1.0 over a wide range of a values, including the second ballooning stable 
domain. Growth rates 7 (a) and frequencies lo v (Jo) obtained using AWECS (solid line with circles) are 
compared with results obtained by Dong et al., 2004 (DCZJ04) [27] (solid line with dots). In addition, 
the MHD frequencies for lst-order aTAE modes (with o;* p i corrections) are plotted (dash-dotted line). 



the earlier results by DCZ99, where max ~ 16 was used. On the other hand, our convergence study 
in Section 7 will show that # max ^ 40 may be necessary. 

Note that the ESITG branch (rjj = 2.5, a < 0.5) is reproduced accurately; both the frequencies 
and the growth rates agree well. ESITG results are not affected by the different A 2 used by ZC02. 

Based on the agreement with results by ZC02 (major discrepancies were explained) , we will assume 
that AWECS correctly computes the properties of finite-/? ESITG (drift- Alfvcn waves), AITG modes 
and KBMs near the first stability boundary. Discrepancies between the initial value code results and 
those obtained with the eigenvalue approach used by DCZ99 and DCZJ04 remain to be understood. 

5.4 High-/3 instabilities 

The kinetic excitation of Alfvenic instabilities near a£j| and in the second ballooning stable domain 
has been studied by Dong et al., 2004 (DCZJ04) using an eigenvalue solver [27], and by Hirose, Zhang 
& Elia, 1994 (HZE94) using a shooting code [23]. In this section, we compare AWECS results with 
results obtained in these two earlier studies. 

We begin with a discussion of a case studied by DCZJ04, shown in Fig. 7, where the following 
parameters were used 

• Physical parameters: q = 1.5, s = 1.0, e„ = 0.2, fc i = 0.212 {k#pf ons = V2fc i = 0.3), = 1.0, 
r) s = 2.5, e = 0. 

• Numerical parameters: N m = 512 x 4 x 7, max = 40, N g = 512, At = 0.01. 

In the domain scanned by DCZJ04, 0.3 < a < 3.1, there is a significant difference in the growth rate 
and critical a values. In fact, DCZJ04 reported that no instability was observed for a > 3.1; which 
does not agree with our findings. Although, we are not able to verify the growth rates (except for a 
convergence study), we can confirm that the instability observed in our results is a physical mode: The 
frequency in the region 2.8 < a < 4 follows closely the frequency of the first-order aTAE obtained with 
the MHD shooting code, which is denoted by ^tae* m Fig- 7 (dash-dotted line). The corresponding 
mode structure is plotted in Fig. 3(f). 
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Figure 8: Benchmark for s = 0.4 over a wide range of a values, including the second MHD ballooning 
stable domain. Growth rates 7 (a) and frequencies lo t (b) obtained using AWECS are compared with 
results obtained by Hirose et ai, 1994 (HZE94) [23]. In (b), the frequencies of the primary unstable 
mode from HZE94 is plotted. In addition, the MHD frequencies for 1st- and 2nd-order aTAE modes 
(with uj^pi corrections) are plotted (dash-dotted lines). 

Here and in the following, w^| Ej denotes the w* P i-corrected MHD frequency of an aTAE. The 
"order" j ' > 1 identifies in which potential well the mode amplitude has its maximum, and p > counts 
the number of zeros the mode structure has in potential well j [7, 28]. The ground state corresponds 
to p = 0. The frequencies w^tae* are calculated by a MHD shooting code solving Eq. (63). In this 
work, only solutions obtained by shooting along the real 6 axis are considered. A more complete study 
would require shooting into the complex plane and the application of phase-integral methods [28], 
since oTAEs are generally coupled to the Alfvcn continuum, either directly or via barrier tunneling. 
This implies that the problem requires an outgoing- wave boundary condition and that, in general, the 
solutions of interest are not square integrable. This may explain why the eigenvalue solver DCZJ04 
yields a different result. 

For a > 4.8, a 2nd-order aTAE, with frequency w r < w^ Et , becomes the dominant instability. 
The corresponding mode structure is plotted in Fig. 3(g). In the MHD limit, this mode is strongly 
damped in the entire a-range scanned, so no shooting result for w^tae* * s available. However, a fcor 
scan at a = 6.0 revealed that w r (fcoi — > 0) is nonzero, which supports our assertion that the mode is 
an aTAE. 

Let us now proceed to the case studied by HZE94, shown in Fig. 8, where the following parameters 
were used 

• Physical parameters: q = 1.2, s = 0.4, e„ = 0.175, fc i = 0.1, = 1.0, t]\ = t] e = 2.0, e = 0. 

• Numerical parameters: JV m = 512 x 4 x 7, (9 max = 60, N g = 1024, At = 0.02. 

It must be noted that the term v\\dg is omitted in the model used by HZE94, so that important physical 
effects associated with the transit resonance are missing. 

Figure 8(a) shows that the qualitative behavior of the growth rate seen by HZE94 is reproduced 
by AWECS, both for KBM and aTAE. The fact that the growth rate in HZE94 are systematically 
larger is most likely related to their neglect of the transit resonance (Landau damping). Note that 
the transition to a 3rd-order mode around a « 2.8 apparent in HZE94's results is not obvious in the 
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Figure 9: Benchmark of AWECS against ATAE for resonant excitation of aTAEs by trapped energetic 
ions. In (a), the ballooning potential Vbaii and the mode structure |<$\l/ s (0)| = y/J\S^ c (6)\ are shown. 
In (b), the instantaneous contributions of individual phase-space markers to the kinetic compression 
term are plotted as a function of the resonance condition K = — uJd)/wb, distinguishing between 
particles with positive (+) and negative (reversed) precessional drift (o). Most non-resonant particles 
are screened out in (b). The snapshot is taken at t = 400 and the growth rates and frequencies are 
average values for t £ [320,400]. The errors represent the standard deviation due to discretization 
noise; the uncertainty of the mean values is smaller than this noise level by a factor v81 = 9. 



AWECS data. For more accurate simulations in this regime, markers must be loaded in a broader range 
of the simulation domain. The size of the simulation domain itself may also need to be increased. 

Near the second MHD ballooning stability boundary, in the range 1.5 < a < 1.6, the mode 
frequency u> T is similar to that of the lst-order aTAE obtained with the MHD shooting code, w^tae* ' 
as can be seen in Fig. 8(b). In the range 1.6 < a < 3, a 2nd-order aTAE is excited. This mode is 
strongly damped for a < 2.8, so shooting results for w^ T ae* are available only for a > 2.8. However, a 
fcorscan at a = 2.3 revealed that w r (fcoi — ► 0) is nonzero, which supports our assertion that the mode 
is an aTAE. 

In summary, the qualitative agreement between AWECS results and those by HZE94, and the agree- 
ment of the frequencies with those obtained for MHD aTAEs suggests that AWECS accurately repro- 
duces essential properties of high-/3 shear- Alfven instabilities. Discrepancies with results by DZCJ04 
remain to be resolved. Further simulations using AWECS indicate that the effects of SBn and fl p are 
small in the a range scanned in the two cases discussed above. As an example, several points obtained 
using the full code are plotted in Fig. 8 (bold circles). 



6 Benchmark 3: Resonant excitation of aTAEs by energetic 
ions 

The excitation of aTAEs through trapped energetic ions was previously demonstrated by Hu & Chen, 
2005 (HC05) with the Sf PIC code ATAE [28]. In this section, we use ATAE results as a benchmark for 
AWECS. ATAE employs a hybrid model consisting of an ideal-MHD core plasma (thermal electrons and 
ions) and a sparse population of trapped energetic ions. The dynamics of the latter is governed by the 
GKE (15). The equations are sufficiently simple to be advanced with a leap-frog scheme, which allows 
to solve the GKE directly, without the need to apply the substitution SG — > Sg, Eq. (24). 

In order to produce results comparable to those of ATAE, several adjustments needed to be carried 
out in AWECS: 

• We omit all terms containing A' s ; i.e., 6S2 and SA 2 and (un A' JiSG), which are neglected in ATAE. 
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• For the GKE, we adopt the ^-dependence of Cl K , fl p and w* s as defined in ATAE, where these 
quantities are all oc 1/B. 

• Let k±p C [ <C 1, take the ideal-MHD limit, 8E\\ = 0, and approximate u)5B\\ = —VLpSifj. In ATAE, 
we have replaced 

n p = -( M o/B 2 )fei-(bxVP c ), 

where P c = P — Pe, by 

ftp = -(po/B 2 )^ ■ (b x VP), 
since the ordering (3e/Pc ~ 0{e) used in Ref. [28] does not apply for the parameters used. 

• We eliminate effects of passing energetic ions. For trapped energetic ions, we retain only the 
ballooning term and the kinetic compression. 

The most important modifications are the neglect of thermal ion FLR effects, the parallel electric 
field and passing energetic ions. While the former two are easily implemented by setting fcoi <C 1 and 
SU e = 0, the latter requires careful adjustments. After the parallel Ampere's law is reduced to 

SE e « (1 - r 0i )<5* c « M*e, 
we are left with the vorticity equation 

fk 2 0i d 2 S^ c =kl [f6&: + 2hh'8%] + Zf ( Wd J «?) E>trap (65) 

^ei 

g5* e + Zet&t? (J 2 cj d ^F ) Eti . ap S^ c . 

The terms on the second line of Eq. (65) are components of the so-called ballooning term. Their 
drift-kinetic limit may be summarized as k^agS^c (not done here). The core component is 

4w»iO Ki = fc„i(a - a E )g = k^acg, 

the contribution of passing energetic ions is &Q;aE,pass<7, an d that of the trapped energetic ions is 
contained in the term ^ Jq lu^luJ Fq^ e t . The quantity cue. pass is calculated using the phase-space 
marker distribution for passing energetic ions at t = 0. Since the velocity space moment (■■■Fq) in 
Eq. (65) involves only trapped ions it needs to be evaluated numerically. Note that 

r 1,2 (f5K)' = <J< - V hall S% - a(g/f)S% 

[cf. Eq. (63)], where the last term, a(g/f)5^ s , cancels with the drift-kinetic limit of the ballooning 
term. 

In the benchmark simulation, the following parameters are used: 

• Physical parameters: q = 2.0, s = 0.5, a = 2.1, e„ c = e„ E = 0.2, v t E — 0.7, fcoE = 0.21, = 1.0, 
rj c = 0.0, r/E = 1.0, e = 0.2, Ze = Me = 1. Some relevant parameters derived from the above 
are: = 5 x 10~ 2 , rf; = 392, fc oi = 0.010606, # = 2.5 x 10" 3 , a E = 2.0, a E ,tra P = 1-22. 

• Numerical parameters: N m = 512 x 4 x 3, max = 60, A^ g = 1024, At = 0.04. Since ATAE uses a 
2nd-order leap-frog scheme the time step adjusted to At = 0.02. 

This corresponds closely to the case shown in Fig. 1 of Ref. [28], with the difference that we use bounce 
angles in the range 6 b € [0.1°, 179.9°]. 

The results are shown in Fig. 9. The mode structure in Fig. 9(a) shows an oTAE with dominant 
(1,0) component. In Fig. 9(b), the instantaneous contributions of individual phase-space markers to the 
kinetic compression term are plotted as a function of the resonance condition K = (lo t — ZUd)/wb, with 
ZJd being the bounce-averaged magnetic drift frequency and Wb the bounce frequency [13]. Resonances 
can be observed at even values of K which is consistent with the fact that oTAE(l,0) is an eigenfunction 
with even parity. The precessional drift resonance, K = 0, is due to particles with ZUd > 0, whereas 
the drift-bounce resonances, K > 2, are due to particles with reversed drift Ud < 0. 

The good quantitative (5^ c , ui r and 7) and qualitative (wave-particle resonance) agreement between 
results obtained with AWECS and ATAE shows that the interaction between trapped energetic ions and 
aTAEs is accurately reproduced by AWECS with the reduced vorticity equation (65). Our preliminary 
studies indicate that the inclusion of thermal ion FLR effects, SE11 7^ and passing energetic ions 
could significantly modify the results. To our knowledge, such a case has not been studied before and 
AWECS will be used to advance into this area. 
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Figure 10: Convergence study for a case from Fig. 6 with a = 0.7 and r\\ = 0; i.e., a KBM without 
ITG-drive. (a) shows convergence with the simulation domain size # ma x, and (b) convergence with the 
number of markers N m = A^jpass x 4 x A p (controlled by iVf i pass ). 
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Figure 11: Convergence study for a case from Fig. 6 with a = 0.5 and r/i = 2.5; i.e., at the instability 
threshold of AITG modes, (a) shows that convergence with the simulation domain size 6! m ax is not 
achieved. In (b), convergence with the number of markers N m = Af^pass x 4 x A p (controlled by 
-^Vf.ipass) is demonstrated for a given domain size max = 40. 
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Figure 12: Convergence study for a case from Fig. 8 with a = 2.3 and ryi = 2.0; i.e., a predominantly 
2nd-ordcr aTAE driven by ITG. (a) shows convergence with the number of periods A p in which markers 
are loaded, and (b) convergence with the number of markers N m = Af^pass x 4 x A p (controlled by 

Af ,ipass) • 
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7 Numerical convergence 



Examples for numerical convergence with the domain size # ma x, the number of markers N m = Nf x 
4 x N p , and the number of marker loading periods N p are shown in Figs. 10, 11 and 12. The error bars 
indicate the level of discretization noise, which decreases with increasing number of markers. Only 
cases with e = and without energetic ions are considered here; i.e., N p = Api pass (passing thermal 
ions only). 

In Fig. 10, a convergence study is shown for a KBM in the case considered by Dong et ai, 2004 
(DCZJ04) [27] and Zhao & Chen, 2002 (ZC02) [25] (cf. Fig. 6). DCZJ04 used max « 23 and ZC02 used 
#max < 50. ..100, while Fig. 10(a) indicates that at least # max ~ 40 is required in AWECS for numerical 
convergence. It is thus possible that the results by Dong et al. [26, 27] are not fully converged with 
respect to the simulation domain size. The AWECS results in Fig. 6 were obtained with N m = 512x4x5 
markers, and the equivalent number of markers used in ZC02's calculations is N m « (315. ..630) x 4 x 5, 
which appears to be sufficient according to Fig. 10(b). 

Near the stability boundary, convergence with # max becomes problematic, as is illustrated in 
Fig. 11(a) for an AITG mode [25] (cf. Fig. 6). This result shows that it may be difficult to ob- 
tain an accurate value for a£ it (and possibly as well) other than through extrapolation from 
converged results away from the stability boundary. Despite the lack of convergence with respect to 
#max, the results still converge well with the number of markers for a given (9 max , as can be seen in 
Fig. 11(b). 

In Fig. 12, a convergence study is shown for higher-order aTAEs (with dominant 2nd-order mode) 
in the 2nd stable domain in a case considered by Hirose et ai, 1994 (HZE94) [23]. Due to the broad 
mode structure, markers must be loaded in at least N p = 7 periods, as was done in Fig. 8 from which 
this case was taken. The number of markers was N m = 512 x 4 x 7, which is also found to be sufficient. 

Given the limitations of the model used, the accuracy of the results shown in Figs. 6-8 may be 
considered to be sufficient to delineate the qualitative features of the modes studied. For this purpose, 
it is usually not necessary to carry out the expensive calculations required for full convergence. 

8 Conclusions and discussions 

A 1-D linear gyrokinetic code called AWECS has been developed to study the kinetic excitation of 
Alfvenic instabilities in a high-/? tokamak plasma. The model equations and the numerical scheme are 
described and the code has been tested carefully. In particular, it is shown that AWECS reproduces 
successfully essential properties of ESITG and AITG modes, KBMs, and aTAEs. Benchmarks against 
results in Refs. [16, 25, 19, 23, 28] are regarded as successful. However, discrepancies persist in 
comparisons with Refs. [26, 27]. 

While the real frequencies calculated by AWECS have also been confirmed by MHD shooting code 
calculations, the quantitative accuracy of the growth rates is more difficult to show, and it is here where 
the main discrepancy with Refs. [26, 27] lies. The code used by Zhao & Chen [25] is most directly 
comparable to AWECS (without energetic ions), and here the growth rates agree well. Numerical con- 
vergence studies in several typical cases indicate that accurate calculations require a larger simulation 
domain than used in Refs. [26, 27]. Furthermore, outgoing boundary conditions are required to accu- 
rately reproduce properties of modes subject to continuum damping, such as aTAEs. Hence, the lack 
of numerical convergence and sensitivity to the boundary conditions may explain the discrepancies, as 
was previously suggested in Ref . [25] . 

Overall, it can be concluded that AWECS is functioning properly and that it can be used for further 
research. Since benchmarking cannot rule out all possible modeling and programming errors, AWECS 
will undergo continuing scrutiny while in operation. Code maintenance and extensions are simplified 
by the modular structure of the code and the application of principles of object-oriented programming. 
Interactive graphical tools were developed using matlab in order to assist the user in data analysis 
and post-processing tasks. 

The demonstration of aTAE excitation in the cases studied by Hirose et ai, 1994 [23] and Dong 
et al., 2004 [27] constitutes the first successful application of this code. Note that, in these earlier 
works, the observed instabilities were not identified as aTAEs, which were discovered more recently 
[7] . Details about the physics of aTAE excitation through wave-particle interactions with thermal and 
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energetic ions will be reported elsewhere. With the use of AWECS, many other interesting problems 
may be addressed in the future, including the physics of BAEs, TAEs and EPMs. 

In this paper, the CHT s-a model equilibrium is adopted in order to be able to carry out com- 
parisons with previous studies, and, thereby, benchmark the code. Currently, work is underway to 
implement the local equilibrium model developed by Miller et al. [29] which will allow us to explore 
the high-/? regime more relevant to experiments. Further extensions are under consideration, such as 
the inclusion of effects due to magnetically trapped electrons. 

Finally, two cautionary notes remain to be added. First, since the CHT s-a model is derived under 
the assumption that B is independent of 9, there is a certain degree of arbitrariness in the way how the 
extension to a model with finite aspect ratio, e > and, thus, variable B(9), is carried out. Second, all 
velocity integrals not arising from the transformation SG — > Sg must be evaluated analytically or using 
phase-space markers distributed uniformly along 9. On the other hand, those integrals arising from the 
transformation SG — > 5g must be evaluated numerically, utilizing the same markers used to calculate 
Sg (nonuniform in 9). If this is not done, inconsistencies arise due to a mixing of ^-independent and 9- 
dependent densities. Preliminary tests where the ballooning term contains B, so that the cancellation 
with the ag/f term arising from the substitution S^ c — * S^ s [cf. Eq. (62) and the note following 
Eq. (65)] is incomplete, have shown that such inconsistencies may modify the ballooning potential to 
such a degree that unstable modes appear in the second MHD ballooning stable domain even without 
wave-particle interactions. In summary, care must be exercised when the CHT s-a model employed in 
this paper is used with finite aspect ratio, e > 0. 
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Erratum 

In the published version of this work [A. Bierwage & L. Chen, Communication of Computational 
Physics 4, 457 (2008)], the following errors were found: 

1. The coefficient in front of SU on the left-hand side of Eq. (2.33) was + To; — Hj) but 
should read 

[l/rl + r 0i - (H u - Z e t&t£)] . (66) 

2. The coefficient in front of 8C e on the left-hand side of Eq. (2.36) was (B 2 + A^A U ) but should 
read 

{& +%AlA u ). (67) 

3. The title of Section 4.2 was "Zero aspect ratio" but should read "Zero inverse aspect ratio." 
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